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Some modifications of the Liischer algorithm, which reduce the autocorrelation time, are proposed and tested. 
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1. HUBBARD MODEL IN THE 
LUSCHER REPRESENTATION 

The problem of dynamical fermions stiU at- 
tracts new interest. Although solved in principle 
still provides the considerable challenge in prac- 
tice. To remind, the main difficulty with simulat- 
ing the theory with dynamical fermions consists 
of the nonlocality of the fermionic determinant 
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and o,> denotes the creation 
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exp(-*M*) ~ det{M), 
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which depends functionally on the gauge field in 
the case of QCD for example. The numerical 
complexity of best known algorithms, namely hy- 
brid Monte Carlo, scale with the volume of the 
system as V^/^ [1]. H owever these algorithms 
suffer from the strong correlations between gener- 
ated configurations. Therefore the Liischer's pro- 
posal of algorithm with numerical complexity V 
is very promising [2, 3]. 

In this report we present some new results of 
applying the Liischer method to the Hubbard 
model. Our aim at present stage is to better un- 
derstand strong correlations, existing also in this 
approach [4], and possibly to offer some remedy 
of this problem. Our findings are rather general 
and can be relevant also for other systems. 

The Hamiltonian of the Hubbard system reads 
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operator of an electron at the lattice site i and 
with the spin cr =t, |. The physical parameters 
entering the model are: K - hopping parameter, 
U - the strength of the effective Coulomb inter- 
action, ji - the chemical potential and the inverse 
temperature /3. The detailed derivation of the 
Liischer representation of the Hubbard model can 
be found elsewhere [5] . The final form of the par- 
tition function reads 
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Where Q = M.'^M./Xmax and Xmax denotes 
the largest eigenvalue of M.'^ M., M. being the 
original fermion matrix entering the Euclidean 
formulation. 
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A{x) is the Hubbard-Stratanovich continuous 
field, and 4'k{x) are the auxiliary bosonic fields 
introduced by Liischer, ak and Pk being the real 
and imaginary parts of the zeros of the polyno- 
mial Pjv(a;) approximating the inverse Pjv(a;) ~ 
1/x. 

The Liischer approximation is uniform in the 
interval (e, 1) and the error falls exponentially 
with the number of bosonic fields. The value of 



1 



e is determined by the width of Q^Q spectrum. 
We have kept the relative error at the level of 
0(10"^). This required 100 - 150 bosonic fields 
in the simulations. 

2. RESULTS AND DISCUSSION 

The runs have been done mainly on the 6^8 lat- 
tice. One MC step consists of the single heat bath 
generation for (j) fields and ten Metropolis updates 
of the Hubbard-Stratonovich field A. Addition- 
ally, in order to check our results, the algorithm 
which performs the exact numerical evaluation of 
determinant has been implemented. 

It has been previously reported, in the QCD 
case, that simulations with Liisher technique suf- 
fer from the long autocorrelation [4] times. The 
source of these correlations was also suggested 
there. It is however important to learn whether 
they can be reduced by the selective refinement of 
the generation of the (j) fields or while an update 
of the A field. Below we discuss some improve- 
ments. 

Tab.l gives the results for the two observables: 
the average density of electrons < ti| > and 
that of pairs of electrons with the opposite spin 
< Ti|(a;)Ti|(a;) > . The simple algorithm described 
above gives large autocorrelation times already 
for K = \ and {7=1. Partly they could be caused 
by the critical slowing down introduced by the (jy^ 
fields, especially those with small (3^. Since the 
model is in fact gaussian in (j) fields, it is natural 
to apply multigrid (MG) methods [6]. Indeed the 
MG generation of (j) fields reduces the autocorre- 
lation times substantially (cf. row 3) . However it 
introduces additional CPU coast which is practi- 
cally prohibitive on larger lattices. The maximal 
decorrelation of the auxiliary fields is achieved by 
the independent generation of the eigenmodes of 
the quadratic form , see row 4. This requires the 
inversion of the N matrices and the results are 
comparable to the MG case with W-cycle. 

Neither simple version nor MG refinement is 
capable to reproduce the exact determinant re- 
sults for {7 = 2 (cf. column 3 of Table 1). 
The system did not thermalize even after ~ 20 
times more thermalization steps then required for 
{7 = 1. This fact can be simply understood. 
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Figure 1. The distribution of the conditioning 
numbers for the M.'^ Ai matrix (solid line) and 
preconditioned matrix M.^ Ai (dashed line) 



When {7 increases the spectrum of matrix M.'^ Ai 
becomes exponentially wider and the bigger num- 
ber of (j) fields in polynomial approximation is 
needed. The constraints imposed by large num- 
ber of (j) fields on one A field becomes more re- 
strictive and the mobility of algorithm rapidly 
decreases. Performing updates of A fields more 
frequently is only a partial solution. 

More promising is the direct preconditioning of 
the bosonic matrix. We define the precondition- 
ing procedure as follows. Let D be the diagonal 
part ofM. = A + D. Then the matrix D''^ M 
is better conditioned as can be seen from Fig.l. 
This effect is also clearly visible in the last row 
of Table 1. Our preconditioning reduces the au- 
tocorrelation time by more than a factor of 8. 
As a consequence wider range of couplings and 
lattice sizes becomes available. For largest lat- 
tice correlation time had large errors due to the 
shorter (relative to the Tint measurement time). 
As the conditioning number of M.'^ Ai is largely 
governed by factor (iU jNt one can reduce the 
conditioning numbers by choosing large enough 
Nt. Our recent simulations performed on the 6^14 
lattice confirm this expectations. However still a 
lot of work remains to be done in order to reach 
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Table 1 





3 = 1,U = 1 


3 = 1,U = 1 


3 = 1,U = 2 


3 = 1,U = 1 




lattice 5^ 


lattice 6^8 


lattice 6^8 


lattice 6^14 


Exact determinant 


0.460(2) 


0.473(2) 


0.462(4) 






0.2197(2) 


0.2203(4) 


0.195(1) 






TZ;^ + — 2i 

' %71Z ^ 


TZ;^ + — 2i 

' %71Z ^ 






Simple program 


0.461(5) 
0.2180(8) 

Tint = 570 


0.471(6) 
0.221(1) 

Tint = 870 
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Simple program 
with MG (W-cycle) 




0.470(4) 
0.2194(6) 

Tint = 160 






Global generation 


0.449(7) 








of gaussian fields 


0.221(1) 

Tint = 170 








Simple program 




0.470(5) 


0.463(7) 


0.475(7) 


with preconditioning 




0.220(1) 


0.195(1) 


0.186(3) 






Tint = 100 


Tint = 200 


Tint = 600 



Results for the density of electrons (upper) and those for the density of pairs (lower). The integrated 
autocorrelation times Tint obtained from both observables are similar. They are quoted in units of the 
single update (sweep) of the A field. 



the physically interesting values of the couplings 
while keeping A4^ A4 reasonably well conditioned 
pU/Nt ~ 0.5 [7]. 

The recent modification of Liischer method [8] 
is very promising. Additional global Metropolis 
step, which makes the approach exact, allows also 
better control of the number of auxiliary fields. 

In conclusion, the Hubbard Model is in the 
Liischer class, i.e. it can be mapped onto a system 
of bosons with local interactions. The half filled 
case has the positive Boltzman factor which ad- 
mits standard Monte Carlo techniques. Precon- 
ditioning of the fermionic matrix strongly reduces 
the large autocorrelation times. Further studies 
are needed to turn this approach into viable al- 
ternative to the existing methods. 
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